function [Q,theta2,params,jointDownScalingOn] = induceStatVAR_threshold_macro_spec(paramsValues,input)

jointDownScalingOn = 0;
if any(paramsValues>1)
    Q = 1D35;
    return
end
% We turn paramsValuesinto a struct 
for i=1:size(input.selectParams,1)
    name = input.selectParams(i,1);
    paramsInput.(name{1}) = paramsValues(i,1);
end

% Accounting for calibrated coefficients, we unfold paramsInput by to get params
params = paramsInput;
% The setdiff implies that we do not overwrite elements in params based on calibrateparams
if ~isempty(input.calibrateParams)
    namesCalibrate = setdiff(fieldnames(input.calibrateParams),input.selectParams);
    for i=1:size(namesCalibrate,1)
        name = namesCalibrate(i);
        params.(name{1}) = input.calibrateParams.(name{1});
    end
end

% We unfold theta2
nx = size(input.varXData,1);
[~,~,hx_1,hx_2,~,gama0,gamaz,gamax,sigmaz] = unfoldTheta2_threshold_macro(input.theta2,nx);
% Scaling hx_1 and hx_2
hx_1 = params.delta1*hx_1;
hx_2 = params.delta2*hx_2;
hx_2(input.switch_hx == 0) = hx_1(input.switch_hx == 0);

% Check if stationarity is possible by only downscaling either hx_1 or hx_2
eigValues_1         = eig(hx_1);
modulus_1           = sqrt(real(eigValues_1).^2 + imag(eigValues_1).^2);
eigValues_2         = eig(hx_2);
modulus_2           = sqrt(real(eigValues_2).^2 + imag(eigValues_2).^2);
if any(modulus_1>=1) || any(modulus_2>=1)
    %disp('Downscaling both hx_1 and hx_2')
    hx_1 = min(params.delta1,params.delta2)*hx_1;
    hx_2 = min(params.delta1,params.delta2)*hx_2;
    jointDownScalingOn = 1;
end

% Re-estimating h0_1, h0_2 and sigma
T        = size(input.econAct,1);
regime_1 = input.econAct >= input.thresholdLevel; %a boom
regime_2 = input.econAct < input.thresholdLevel;  %a recression
if input.h0RegimeOn == 1
    h0_1 = (eye(nx)-hx_1)*sum(input.factors(1:nx,:).*repmat(regime_1',nx,1),2)/sum(regime_1);
    h0_2 = (eye(nx)-hx_2)*sum(input.factors(1:nx,:).*repmat(regime_2',nx,1),2)/sum(regime_2);
else
    h0_1 = mean(input.factors(1:nx,:),2)-hx_1*sum(input.factors(1:nx,:).*repmat(regime_1',nx,1),2)/T ...
                                        -hx_2*sum(input.factors(1:nx,:).*repmat(regime_2',nx,1),2)/T; 
    h0_2 = h0_1;
end
h0_2(input.switch_h0 == 0) = h0_1(input.switch_h0 == 0);
Wxhat               = input.factors(1:nx,2:T) - (repmat(h0_1,1,T-1) + hx_1*input.factors(1:nx,1:T-1)).*repmat(regime_1(1:T-1)',nx,1) ...
                                              - (repmat(h0_2,1,T-1) + hx_2*input.factors(1:nx,1:T-1)).*repmat(regime_2(1:T-1)',nx,1);
VarWxHat             = zeros(nx,nx);
for t=1:T-1
    omegaHat = input.VarU(1:nx,1:nx,t+1)...
        + hx_1*input.VarU(1:nx,1:nx,t)*regime_1(t,1)*hx_1' ...
        + hx_2*input.VarU(1:nx,1:nx,t)*regime_2(t,1)*hx_2' ...
        - input.Cov10U(1:nx,1:nx,t+1)*regime_1(t,1)*hx_1' - hx_1*input.Cov01U(1:nx,1:nx,t+1)*regime_1(t,1) ...
        - input.Cov10U(1:nx,1:nx,t+1)*regime_2(t,1)*hx_2' - hx_2*input.Cov01U(1:nx,1:nx,t+1)*regime_2(t,1);    
    VarWxHat = VarWxHat + 1/(T-1-2*(nx+1))*Wxhat(:,t)*Wxhat(:,t)' - 1/(T-1)*omegaHat;
end
[sigma,checkWx] = chol(VarWxHat,'lower');
if checkWx ~= 0
    %sigma = sigmaOld;
    VarWxHat = zeros(nx,nx);
    for t=1:T-1
        VarWxHat = VarWxHat + 1/(T-1-2*(nx+1))*Wxhat(:,t)*Wxhat(:,t)';
    end
    sigma = chol(VarWxHat,'lower');
end


% We have a stationary model
if input.hxRegimeOn == 0 && input.h0RegimeOn == 0
    % Unconditional variance is available in closed-form
    h0_1  = (eye(nx)-hx_1)*mean(input.factors(1:nx,:),2);
    T      = size(input.factors,2);
    residualsTmp = zeros(nx,T);
    for t=2:T
        residualsTmp(:,t) =  input.factors(1:nx,t) - (h0_1 + hx_1*input.factors(1:nx,t-1));
    end
    varU     = sum(input.VarU(1:nx,1:nx,1:T-1),3);
    cov10U   = sum(input.Cov10U(1:nx,1:nx,2:T),3);
    cov01U   = sum(input.Cov01U(1:nx,1:nx,2:T),3);
    omegaTmp = 1/(T-1-nx-1)*(residualsTmp*residualsTmp') + 1/(T-1)*(-varU - hx_1*varU*hx_1' + cov10U*hx_1' + hx_1*cov01U);
    [sigmaTmp,check] = chol(omegaTmp,'lower');
    deScaling     = 1;
    while check ~= 0
        deScaling = deScaling*0.9;
        omegaTmp = 1/(T-1-nx-1)*(residualsTmp*residualsTmp') + 1/(T-1)*deScaling*(-varU - hx_1*varU*hx_1' + cov10U*hx_1' + hx_1*cov01U);
        [sigmaTmp,check] = chol(omegaTmp,'lower');
    end
    varX     = reshape((eye(nx^2)-kron3(hx_1,hx_1))\reshape(sigmaTmp*sigmaTmp',nx^2,1),nx,nx);
else
    % The unconditional variance is computed by simulation
    numSim         = input.numSim;
    thresholdLevel = input.thresholdLevel;
    rng(1,'twister')
    shocks         = randn(nx,numSim);
    shocksz        = randn(1,numSim);
    xSim           = zeros(nx,numSim);
    zSim           = zeros(1,numSim);
    for s=2:numSim
        if zSim(1,s-1) >= thresholdLevel
            xSim(:,s) = h0_1 + hx_1*xSim(:,s-1) + sigma*shocks(:,s);
        elseif zSim(1,s-1) < thresholdLevel
            xSim(:,s) = h0_2 + hx_2*xSim(:,s-1) + sigma*shocks(:,s);
        end
        zSim(1,s) = gama0 + gamaz*zSim(1,s-1) + gamax'*xSim(:,s-1)+ sigmaz*shocksz(1,s);
    end
    varX = zeros(nx,nx);
    for i=1:nx
        varX(i,i) = var(xSim(i,input.burnin+1:numSim),1,2);
    end
end
% Matching all the variances
Q = 0;
for i=1:nx
    Q = Q + ((varX(i,i)-input.varXData(i,i))/input.varXData(i,i))^2;
end
Q = sqrt(Q/nx)*100;


% We save the potentially scaled elements in theta2
theta2 = setTheta2_threshold_macro(h0_1,h0_2,hx_1,hx_2,sigma,gama0,gamaz,gamax,sigmaz);

end

